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Algorithm in Asynchronous MIMO OFDM Signals 
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Abstract — Spectrum sensing (SS) in cognitive radio (CR) sys- 
tems is of paramount importance to approacli tlie capacity limits 
for tlie Secondary Users (SU), wliile ensuring tlie undisturbed 
transmission of Primary Users (PU). In tliis paper, we formulate 
a cognitive radio (CR)systems spectrum sensing (SS) problem 
in which Secondary Users (SU), with multiple receive antennae, 
sense a channel shared among multiple asynchronous Primary 
Users (PU) transmitting Multiple Input Multiple Output (MIMO) 
Orthogonal Frequency Division Multiplexing (OFDM) signals. 
The method we propose to estimate the opportunities available to 
the SUs combines advances in array processing and compressed 
channel sensing, and leverages on both the so called "shrinkage 
method" as well as on an over-complete basis expansion of 
the PUs interference covariance matrix to detect the occupied 
and idle angles of arrivals and subcarriers. The covariance 
"shrinkage" step and the sparse modeling step that follows, 
allow to resolve ambiguities that arise when the observations are 
scarce, reducing the sensing cost for the SU, thereby increasing 
its spectrum exploitation capabilities compared to competing 
sensing methods. Simulations corroborate that exploiting the 
sparse representation of the covariance matrix in CR sensing 
resolves the spatial and frequency spectrum of the sources. 



I. Introduction 

Generally, spectrum-sensing methods include matched filter 
detection HI, |2^|, likelihood ratio test (LRT) 11], energy 
detection Q-lSl, and cyclostationary feature detection Q- 
ifm. each of which has different requirements and advan- 
tages/disadvantages. The non-coherent energy detector has 
been shown to be optimal if the cognitive devices have 
no a priori information about the features of the primary 
signals except local noise statistics. In addition, it obviates the 
need for synchronization with unknown transmitted signals. 
Matched filter based detection requires a priori knowledge 
of the primary user, e.g., modulation type, preambles, pilots, 
pulse shaping, and synchronization of timing and carrier If the 
modulation schemes of the primary signals are known, then 
the cyclostationary feature detector can differentiate primary 
signals from the local noise by exploiting certain periodicity 
exhibited by the mean and autocorrelation of the correspond- 
ing modulated signals. The cyclostationary detection needs to 
know the cyclic frequencies of the primary signal, which may 
not be available to the secondary users in practice. 

These methods detect the presence of the PU within a 
band and they are presumed to be combined with stochastic 
control algorithms that decide strategically what bands the CR 
receiver should examine next lfT2l . Finding spectrum holes in a 
wideband signal is, instead, the objective of wideband SS and 



the focus of this paper The exemplary scenario we envision 
for the CR is that of a femto-cell access point 1 14|, in the role 
of the SU, searching for spectrum opportunities in a dedicated 
band with base stations acting as PUs and transmitting MIMO- 
OFDM signals. OFDM is the modulation of choice for most 
of the emerging broadband wireless communication physical 
layer standards. In [131 the authors argue that OFDM is the 
best physical layer candidate for a CR system since it allows to 
modulate signals so as to fit into discontinuous and arbitrarily 
wide spectrum segments. OFDM is also optimal from the 
viewpoint of approaching the Shannon channel capacity in 
a wideband channel with frequency selectivity and colored 
noise. Finding the spectrum opportunities in an OFDM system 
is equivalent to detecting the spectrum holes due to unoccupied 
subcarriers. A similar problem has been formulated in (ITSl . 
|16|. The basic difference in our model is that we do not 
assume any form of OFDM symbol or frame synchronization 
among PUs and between PUs and SUs. This reflects the wide 
practice, for example, in today UTE systems and also allows 
to tackle more general applications. 

Usually, the PUs do not cooperate with the SU and do not 
transmit specific synchronization signals for the purpose of 
synchronization at the SU. Synchronization sequences from 
PUs are going to be sent in each frame to allow their 
mobile users to synchronize with the frame period. We assume 
that the CR has lined up to the strongest PU signal frame 
period and focus instead on the rapid estimation of spectrum 
holes, via second order methods by estimating the interference 
covariance matrix. The structure of the covariance matrix 
is dictated by the cyclic prefix in OFDM symbols, antenna 
array manifold, occupied subcarriers and channel parameters. 
In this paper, we focus on the rapid estimation of spectrum 
holes, via second order methods by estimating the interference 
covariance matrix. The motivation behind using a second 
order method is that it is non-coherent, does not require 
synchronization, and does not need knowledge of the PUs 
modulation. 

Previous papers that focused on wideband SS for CRs, 
used detectors that leverage the structure of the second order 
statistics of the PU signal 1 17|-|21 1. Compared to these papers, 
our MIMO-OFDM sensor approximates a Generalized Likeli- 
hood Ratio Test (GLRT) using an estimation algorithm for the 
covariance that we call Shrink and Match (S&M). The S&M 
algorithm approximate the maximization of the likelihood 
function for Gaussian PUs, with respect to the parameters of 
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the PUs signal covariance, by alternating between a step of 
the shrinkage algorithm f24\ and a step of the Orthogonal 
Matching Pursuit (OMP) algorithm 1,31] on an appropriately 
defined linear sparse model for the interference covariance. 
Under appropriate conditions, the shrinkage method converges 
to the ML estimate of the Gaussian PUs covariance [25 1, [26] 
with a relatively low cost iteration, and the OMP algorithm is 
used to denoise the estimate obtained at each step, leveraging 
the sparsity of the model. 

A widely explored trade-off in SS is between the time used 
for sensing and that used to exploit the channel (22^. Thus, 
CR SS must work with very short data records. Harnessing 
the benefits of its two steps, the numerical results shown the 
S&M method features excellent performance in this regime. 

Notation: The set of real, complex and integer numbers 
numbers by R, C and Z, respectively. We denote sets by 
calligraphic symbols, where the intersection and the union of 
two sets A and B are written Ar\B and AiJB, respectively. 
The operator |^| on a discrete set takes the cardinality 
(measure) of the set and A"^ denotes the complement of A, 
where the universal set should be evident from the context. 
We denote vectors and matrices by boldface lower-case and 
boldface upper-case symbols. The transpose, conjugate, Her- 
mitian (conjugate) transpose, inverse and pseudo inverse of 
a matrix X are denoted by X^, X*, X^, X^^ and X^, 
respectively. |X| and tr(X) denote the determinant and trace of 
matrix X, respectively. In this work, the vectorization operator 
for a matrix is denoted by vec(X). [X]a ([x]a) is the ath 
column (element) of X (x). Similarly, [X]^ ([x]^) is defined 
as the collection of columns [X]a (entries \x\a) where a ^ A. 
The Frobenius norm of a matrix is denoted by l|X||i?. The 
conventional £2 -norm is written as ||x||2 and ||x||o is the 
number of non-zero entries of the vector x. ® denotes the 
Kronecker matrix product and denotes the Hadamard matrix 
product. The operator E{ } denotes the expectation operator 
and a circular symmetric complex Gaussian random vector 
X G with mean //, G C" and covariance matrix S G C"^" 
is denoted as x ^ CAf{fi, S). 1^ is a w x w matrix of ones and 
I„j is the identity matrix of size m. Omxn denotes an m x n 
matrix of zeros. 

II. System Model 

In our setting, the set of independent asynchronous multiple 
antenna PUs is denoted by I where / = \I\ is the maximum 
number of active sources. The number of available subcarriers 
is denoted by N and Lp represents the cyclic prefix length. 
The length of one OFDM symbol is denoted by M ^ N + Lp. 
The parameter T — 1/W is the sampling period and W is the 
bandwidth of the system. Each PU uses a set of subcarriers 
Ci for transmission where Ci C C — {0, 1, . . . , — 1}. The 
transmitted stream of discrete samples from the ith PU can be 
described as 

Ui [m] = ^ Ui.fe [m - kM] , ( 1 ) 

fcez 

where to G Z and Ui fe[TO] is the k\h OFDM symbol in 
the transmitted data stream from the ith source and can 
be represented in matrix form as Ui ^fm] = Ai ^f/^ with 



TO G M = [0,ilf). The Nt X N matrix Ai^k contains the 
transmitted data symbols in the frequency domain. The vector 
a.i^k,c — [A-i,k]c+i G C^rxi j-jjg modulated data symbols 
transmitted on subcarrier c G Ci in the fcth OFDM symbol 
or fc c — if c ^ Ci. In this work, we assume that the 
vector ai c for c G Ci is randomly distributed zero mean 
with covariance matrix 

E{a,,fc,caffc_J = ^^^I- 

The vector f,'„ = [F],„+i for to G where the A^ x M 
matrix F is constructed by appending the last Lp columns of 
the N X N IFFT matrix at the beginning. 

As mentioned before, we consider a MIMO channel where 
the number of transmit and receive antennas are A^t and 
Nn, respectively. The typical propagation channel in wireless 
systems is assumed to be a multipath channel with at most L 
dominant propagation paths from scatterers in the far field. 
We assume that the antenna spacing is small so that the 
narrowband array manifold approximation holds. In this case, 
a general representation for the discrete time MIMO channel 
between source i and the SU receiver is as follows 

L 

Hj(m, n) = ^ h^^i *(to, eix)g.,{nT - Ti^g) , (2) 
£=1 

where hi^i is the channel fading coefficient with unknown 
variance. The parameter Ti^e is the delay of ith path of 
source i and the matrix ^^(m, 0i_^) is an Nr x Nt MIMO 
channel matrix parameterized by which also depends on 
TO because of Doppler effects and carrier offsets. The function 
gi{-) is the cascade of receive and transmit filter of source i. 

A. Signal Model at the Receiver 

The received signal sampled with the rate of T is: 

y[m] = X! X! Hi(TO,n)ui[TO ~n~ti]+ w[m] , (3) 

i£X n 

where w[to] G C^'^ is a zero-mean Additive White Gaussian 
Noise (AWGN) which is both spatially and temporally white 
and independent of the sources with w[m] ^ CA/'(0, ct^Iat^ ). 
The unknown parameter ti E M. models the misalignment be- 
tween the received OFDM symbols from the ith asynchronous 
PU and the first observation window (starting at m = 0) at 
the SU (See Fig. [B- 

In order to simplify the receiver model, following the 
approach in ||29l , ||30L we discretize the parameters in our 
channel model with a certain resolution. We quantize the delay 
and the parameter vector Oi i and denote the quantized values 
by 

Q{n,i) = q^,tT , ^ Qe(6'^,£) e A& . (4) 

In the following derivations, we ignore the model mismatch 
error The parameter T is the time resolution of the time 
quantization grid and the integer value qi p G Q = [0, Lp) is 
the index of the discrete delay. A@ represents the finite set of 
quantized parameter vectors, whose cardinality is n® = \A@\. 
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Fig. 1 . Contribution of different multipatli components of received OFDM signals from asynchronous PUs to the observed sequence at the SU. 



The function Qe( ) is the quantizer associated with the 
parameter vector Oi^i which can be expHcitly described as 



Qe{Oi,£) = argmin WO^^e - e\\r^ 
eeAe, 



The discrete-time received signal is hereinafter modeled as: 

L 

y[m] = X! X! X! '^tc *("^' 0t.e)ui,k[iT^ - - q,,i - U] 



w[m] , 



(6) 



where the channel coefficients are independent for different i, 
I and fc: 



III. Sparse Covariance Matrix Representation 

In this section, we present a sparse model for the PUs 
MIMO-OFDM interference covariance matrix for the asyn- 
chronous model in (|6]l. In order to accurately estimate the 
covariance matrix, we require multiple {K) independent and 
identically distributed (i.i.d.) samples of the received signal 
with the same covariance structure, which means the param- 
eters {(Ti^f, (7i.f , "'^ remain the same in the time 
required to collect K samples. To ensure that these observation 
samples are i.i.d., the SU should use a sample shift equal 
to 2Af between consecutive collected sequences to guarantee 
that firstly all collected data sequences share the same set of 
parameters {tijf^i and secondly they are uncorrected and 
independent. Thus, the rth NrM x 1 observation sample in 
the SU is expressed as 



y[r2M] 
y[r2M + 1] 

y[r2M + M - 1\ 



(7) 



The following Lemma summarizes the assumptions and the 
mathematical sparse model for the covariance matrix which 
depends on the unknown channel parameters and the unknown 
set of occupied subcarriers. 



Lemma 1: Assume that Eja^jj ^a-i ^ cl 



NT\C^\ 



I for 



c e 



Ci and the transmitted symbols are zero mean and independent 
for different i, k and c. The received signal is zero-mean and 



the sparse representation of its covariance matrix is expressed 
as 



(5) vec(£(^,a2 )) ^ vec(E{y,,yf }) = + al^tciluNn) , 

(8) 

where cr is an MNn@ x 1 sparse coefficient vector with non- 
negative entries. The (MNn)^ x MNn@ matrix M is the 
over-complete dictionary and its columns are constructed with 
vec{Il° g^^ + Ul g J for V e M, 9 e A& and c e C. The 
parameter v capture the unknown delays and misalignments, 
c denotes the subcarriers index and 6 models all the possible 
remaining unknown channel parameters. The matrix YV^ g ^ 
for j — 0,1 is defined as 

^ (J^iv) ® I)r({*(m, 0)}^to')(f,ff I) 

xT{{^''{m,e)}'^=^){{3\v)f®l) , (9) 

[F-'^Jc+i. J-'(w) for j — 0,1 is the shift matrix 
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where fc = 
defined as 



OvxM-v 



Om- 



,3\v) 



Otix Af-i 







AI—vxv 



The operator 7~({A„j}^~q^) creates a block diagonal matrix 
with blocks Ao, . . . , Am-i that appear in its input argument. 
Proof: See Appendix A. ■ 
It is worth mentioning that the sparse model in dHJ for the 
covariance matrix satisfying the constraint cr > 0, guarantees 
that I](cr,cr^y) is Hermitian positive semi-definite. Next, we 
present our proposed method for estimating the structured 
covariance matrix in dHJ. 

IV. The Shrink AND Match Algorithm 

In the CR application, one important objective of the SU 
is to decrease the sensing time in order to increase the time 
to exploit the empty subcarriers for transmission. As a result, 
we focus on the cases where the number of observations K is 
small and the sample covariance matrix S = X]r^=o^ YrY^ 
cannot be considered a good estimate of the true covariance 
matrix. Our proposed method relies on finding first a sub- 
optimal solution that approximates the ML estimate for the 
case of Gaussian PUs symbols and then on matching it to the 
sparse model in Lemma 1 . More specifically, the dxK matrix 
Y = [yo, • • • , Yk-i] is the matrix of collected zero-mean i.i.d. 
samples with the same covariance structure where, in general, 
d denotes the dimension of the problem. Even if the PUs are 



4 



not Gaussian sources, and the ML receiver would suggest to 
perform a complex joint detection, can be approximated 
to be complex Gaussian distributed, as a maximum entropy 
approximation of its distribution. The likelihood function for 
the unknown parameters cr and cr^, given the observation Y 
is expressed as 

exp [^u{YY"^-\cT,al))] 



C{cr,ai\Y) = 



(10) 



7r^'i|S(<T,a2)|i^ 
The sparsity regularized log-likelihood function of (fTOt is 

a,dl=aTgmmlog\-E{(T,al)\+ti\i:-\<T,al) s] + k\\(t\\ 

S.t. crf„>0, (T>0 (11) 

where ||<t||o is the imposed sparsity constraint and k is some 
regularization parameter The ML estimate of the covariance 
matrix is given by 

In ( fTTI ). we have relaxed the original problem by imposing 
the sparsity constraint. However, even solving this relaxed 
problem, which is non-convex, is in general complicated. In 
the following, in order to simplify the optimization, we assume 
that the noise variance cr^, has been estimated and is known. 

Under the small sample size constraint K < d, we ap- 
proximate the solution of (fTTT i by first obtaining an accurate 
estimate of the covariance matrix which is invertible and then 
projecting the estimate of the covariance matrix in the sparse 
model in We call this method the Shrink & Match (S&M) 
algorithm. The two steps are explained next. 

1) Shrinkage Step: For estimating the covariance matrix, 
we use the method proposed in ll24l . In ll24l . the covariance 
estimation is based on shrinkage regularized fixed point iter- 
ations for a high dimensional setting, where the fixed point 
iterations converge to the ML estimate. However, in ll24l . the 
structure of the covariance matrix has been ignored. 

This shrinkage method works with trace-normalized covari- 
ance matrices (tr(E) = d) and in this work, the final estimate 
must be scaled by the estimated value of the trace. Thus, we 
need to estimate the trace of the covariance matrix which is 
obtained as 
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The shrinkage coefficient is estimated as ||24| 

d2 - i tr(RR^) 

d 



so that it is trace-normalized and return to the update step (fT4l) 
until the stopping criterion is met. 

2) Matching Step: In order to impose the structure of 
the covariance matrix, after the convergence of the iterative 
process, we first scale the estimated covariance matrix (as 
S ^ ^Tr(S)Sj_|-i ) to compensate for the trace-normalized 
assumption, and then fit the covariance matrix with the model 
in (O by minimizing the following cost 



d^ - Kd- K 




K 



(K-l) 



tr(RR^) 



(12) 



K <d 

K >d 
(13) 



2 is the trace-normalized 



where R ^ # Et"o' Y.yf /l|y 
sample covariance matrix. The iterative steps update the co 
variance matrix estimate at each iteration as follows 

A'-l 



■'3 + 1 



(1 



H 



y?S^- 



■7l 



(14) 



a = argmin 

cr 

(T > 



vec(S - all) - Mo- 



+ A|k| 



s.t. 



(16) 



At this step, we normaUze the estimated covariance matrix as 

d S ,, 



^3 + 1 



tr(S,+i) 



(15) 



Remark 1: The S&M algorithm solves a problem similar 
to the asymptotic ML (AML) estimator for a structured 
(Hermitian Toeplitz) covariance matrix posed in [32], however 
we have two main differences. In ||32l , the authors use the 
inverse of the sample covariance matrix to define a weighted 
£2-norm that they use to estimate the parameters in their 
model, an approach that is close to optimum in the asymptotic 
regime of K ^ d. Experimentally, we have observed that 
the choice of the shrinkage method to estimate the covariance 
matrix, followed by our sparse denoising step using the £2- 
norm instead of the weighted norm, work better when K < d. 

The problem in (fTSI l can be solved by using greedy methods 
and in particular the non-negative OMP algorithm ll33l in order 
to satisfy the non-negativity constraint on the variables. It can 
also be relaxed by imposing the sparsity regularized constraint 
II • II 1, and in this case, can be solved by convex programming as 
a linear program. Algorithm [T] summarizes the steps required 
for this suboptimal method to approximate the solution of (fTIT i. 
Algorithm |2] summarizes the steps of the Non-negative OMP 
method to solve ( fTSI l. 

Algorithm 1 S&M Algorithm 
Require: Y, M and ct^ 

1: Initialize: j = and Eq — Id- Calculate tr(S) and 7. 

2: Repeat 

3: Calculate Sj+i from ( fT4] i. 



4: Normalization: S 



d S 



"tr(S,+i) 

5: set j = j + 1 

6: Until Stopping criterion: ||I]j— 



l^i-illl 



7: Solve (fT6] l using Algorithm |2] to find a. 



V. The Separable S&M Algorithm 

While the channel model used so far is valid for various 
array configurations and models for time variations, in this 
section we focus on the more familiar case of uniform linear 
arrays (ULA) for the receive and transmit antennas and a 
single Doppler per path. The reason to focus on this case is 
because the covariance has a separable structure that allows to 
introduce a greatly simplified version of the S&M Algorithm. 

The discrete time channel between the SU receiver and the 
ith PU can be described [|27l. Il28l bv 

*(to, e,j) = e^2-^».^'"^e,(/3,,,)ef (a,,,) , (18) 
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Algorithm 2 Non-negative OMP Algorithm 

Require: obtain M = [mi, m2, . . . , m^)], where D 
MNn@ and b = vec(S - a^I) 
1: Initialize: ctq = 0, 5o = 0, ro = b and j = 1. 
2: Repeat 

3: compute 

^2 



(max{mjrj_i ,0})^ 
l|m,||i 



e(d) = l|r,-i||^ 

for 1 < d < D. 

d* = argminrfg5o_^ e(d) 

Update Support: Sj = 5j i U {d*}. 

M, = [M]5, 

Find the Non-negative Solution: 

(T* = argmin ||M,cr - h\\l 



s.t. cr > 



8: Update Solution: 







(17) 



9: Update Residual: = b Mcrj. 

10: set j = j + 1 

11: Until Stopping criterion: ||<Tj — crj_i < ToMp||cj-i | 



where 6ij 



[i>i.i,l3i,e,ai,i]. The vector 64(0^,^) 



] is the steering vector asso- 
ciated with the angle of departure (AoD) and er(/3i.^) = 

]^ is the steering vector asso- 



, ,e 



j27T{NR-l)f:i,^e]T 



[l,e. 

ciated with the angle of arrival (AoA), where the parameters 
ai^i and model the angles of departure and arrival of 
the £th propagation path between the SU and the ith PU, 
respectively. In this model, the frequency offset is incorporated 
into the Doppler spread of the channel. Thus, the parameter 
"04/ — fi + '^i,£ models both carrier frequency offset denoted 
by fi and Doppler effects of the £th path of the ith PU denoted 
by uji^i with ■(/'max — maxi^£ ipij 

The quantized counterpart of Oi^^ is 



(19) 



where the parameter A/ = 1/NT denotes the subcarrier fre- 
quency. The frequency offset is modelled as Pi^i/Pj to capture 
the fractions of the subcarrier frequency (with the resolution 
of A//P/) where the integer value pij e V ^ [pi,pi + P) 
where pi is the smallest value in the set and \V\ — P. 1/A 
and 1/B are the angular resolution of the angle of departure 
and arrival, respectively. Moreover, the integer parameters 
a,^i e ^ = [0, A) and b,^e e S = [0, B). 

Corollary 1: For the quantized discrete time channel model 



e'^^'^er{hj,lB)e^{a,j/A) , 



the (MNji)'^ X MNPB over-complete dictionary M is con- 
structed as 

mn(v,p,c,b) = vec((*p,, T(«)) ® (e,(6/S)ef , 

(20) 

where r]{v,p, c, b) ^ vPNB + (p- pi)NB + cB + b+l. In 
this model, v G A4, the parameter p models the Doppler, 
c € C and b £ B. The matrix T{v) denotes a masking matrix 
and is defined as 



T{v) ^ 



(21) 



The matrix ^p^c — ^{p)^c^^ {p) represents the depen- 
dency of the dictionary components on subcarrier c and 
discretized Doppler p where A{p) = diag{e''^^ '''^ }^Iq is 
a diagonal matrix containing the Doppler coefficients. The 
vector cr in the sparse model dHJ is a MNPB x 1 sparse 
coefficient vector with non-negative entries. 

Proof: See Appendix B. ■ 
In the following, we explicitly explain how the S&M 
algorithm can be simplified in this scenario. In the special 
case where Nji > IL, it is more convenient to divide 
the problem into two separate smaller and easier to solve 
covariance estimation problems, in the spatial and temporal 
domains reducing the memory required to store the dictionary 
elements, the complexity of the algorithm, decreasing its 
runtime. This is possible because of the far field narrowband 
approximation in the channel model, which allows to approx- 
imately separate the effects of path delays and of the different 
arrival times at each array element. We call this method. 
Separable Shrink and Match (SS&M) algorithm. In the first 
step, the SU estimates the AoAs, exploiting the existing 
structure in the spatial covariance matrix that depends on the 
parameters {^i.^j^Z^'' Then using the AoAs information, 
the SU, for each angle, filters spatially the temporal streams 
to produce observation free from interference from the other 
active directions. Then, it uses the temporal sparse covariance 
matrix representation to recover the occupied subcaniers. The 
mathematical formulation of these two steps is given in the 
following. 

Step 1: To estimate the spatial covariance matrix S5 = 
E{y[m]y[m]^}, the SU samples the received signal at the 
antenna array at m = kM and collect K i.i.d. observations 
Ys ^ [y[0],y[M], . . .,y[{K - 1)A/]] with the same spatial 
covariance matrix. The sparse representation of the spatial 
covariance matrix can be expressed as 

vec(Ss) = Mso-s + cr^,vec(lAr^) , 

where erg is a i? x 1 sparse coefficient vector with non-negative 
entries. The N'^ x B spatial dictionary matrix is constructed 
as [Ms]b+i = vec(e^(VS)ef for b £ B. 

Step 2: Similarly, to estimate the temporal covariance matrix 
St, the receiver collects K observation blocks Yxir) — 
[y[2rAf],...,y[2rM + M-l]], r £ {0, 1, . . . , i^T - 1} of size 
Nji X M where 2M is the amount of sample shift between 
consecutive collected sequences. The temporal covariance ma- 
trix of the spatially filtered observations Yrir) = ^^Yrir) 
is expressed as St./ = Ejy^jy* ;} where y,.^; is the kh row of 
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YT{r). $ is the subspace of AoAs and denotes its pseudo 
inverse. The sparse representation of the temporal covariance 
matrix is 

vec(ST,i) Mt(tt,i + fJ^ ||2vec(lA/) , 

where ctt.i is a MNP x 1 sparse coefficient vector with non- 
negative entries. The Af ^ x MNP temporal dictionary matrix 
is constructed as [M.T]f_i(v,p,c) — vec(^'p,c where 
fi{v,p, c) = vPN +ip- pi)N + C+1. 

For the derivations of the sparse representation of the spatial 
and temporal covariance matrices see Appendices C and D, 
respectively. Algorithm |3] summarizes the two steps of the 
SS&M algorithm. 



Algorithm 3 SS&M Algorithm 

Require: Y5, Ms, {YT(r)}lV' Mr and al 
1: Apply Algorithm [T] with M5 and Y5 and find the non- 
negative sparse vector Bs- 
2: Set of detected AoAs: A4, = {b : [Bsh+i > 0, 6 e B}. 
3: Subspace of AoAs: * = {er{b/ B)}b,£A^- 
4: Compute the pseudo inverse of $ 
5: Yrir) = $^YT(r) for r = 0, . . . , ii' - 1. 
6: for / = 1, . . . , 1^0! do 

7: Apply Algorithm[T]with Mt and Y/ = [y^;, . . . , y^_i J, 

and find the sparse non-negative vector Bt,i 
8: £1^ {c: [5^T,/]^(i,,p,c) > 0,^; e M,p eV,cEC} 
9: end for 

10: Set of occupied subcarriers: £ = Ul=i 



VI. Numerical Experiments 

In this section, we examine the performance of our proposed 
algorithm numerically. We consider a MIMO OFDM system 
with 1 = 4 uncorrelated asynchronous sources where Nt = 2 
is the number of transmit antennas. The other channel and 
system parameters are considered to be = 64, Lp = 8 and 
L = 2. We model the first arrival of the OFDM symbols from 
each user as a uniform random variable in Ai and the delays 
are distributed uniformly within the cyclic prefix duration. The 
Doppler pi,i of each path is generated as a uniform discrete 
random variable in the set V with P = 3. The AoAs and AoDs 
are continuous values in [0, 2tt] where the angular difference 
between the AoAs is greater than 10°. In the simulations, we 
have set the number of grid points A = B = 180. Thus, 
the angle resolution of 1 degree. We generate uncorrelated 
Rayleigh fading coefficients, h'j^ ^ ~ CA/'(0, 1). The same 
transmission power equal to 1 is considered for all the sources. 
Throughout this entire section, the signal to noise ratio (SNR) 
is defined as — lOlogj^Qcr^ where the noise is zero-mean 
AWGN with variance equal to cr^. In Algorithm [T] the value 

of Train IS Set tO be 10 

In order to demonstrate the capability of SS&M to make 
efficient use of the empty subcarriers without causing harm- 
ful interference to the PUs, we define two metrics. In the 
following, probability of false alarm and missed detection 
in the cth subcarrier are denoted by Pfa(c) and Pmd(c), 
respectively. The opportunistic spectral utilization of subcarrier 



c is measured with 1 — Pfa(c). As a result, the aggregate 
opportunistic throughput of the SU can be described as 

Af-l 

Rc{1-Pfa{c)), 

where Rc denotes the throughput achievable over the cth 
subcarrier if used by the SU. Assuming that Rc is the same 
for all c G C, the expected aggregate opportunistic throughput 
of the SU can be measured by 

1 

c=0 

Assuming that all the PUs are equally important, the aggregate 
interference to PUs can be expressed as 

N-l 

where Cc denotes the cost incurred by causing interference 
with a PU in the cth subcarrier. Assuming that Cc is equal 
for c G C, the average aggregate interference to PUs can be 
equivalently measured by 

N-l 

Pi -Y ^md(c)- 

c=0 

Ideally, we want to have pr close to 1 while at the same 
time guarantee small values for pi below some threshold e. 
In Fig. 121 pt and pi have been plotted versus SNR when 
Nji = 12 for K = 20 and K = 30 for the same simulation 
runs. The curves are obtained by computing Pmd and Pfa for 
each subcarrier c numerically in 100 independent simulation 
runs. This figure illustrates that SS&M exploits the empty 
subcarriers very efficiently (high aggregate throughput close 
to 1) while causing small interference to the primary users 
(relatively small values for pj) even for small sample sizes. 
In addition, we can observe that increasing the SNR or K 
does not increase pT dramatically. However, increasing K or 
SNR decreases pj. This trade-off is a critical aspect of this 
algorithm that one should consider to choose K based on the 
available SNR and the existing constraint on pi < e. 

To the best of our knowledge, in the literature, there is not 
any method that can be fairly compared with S&M in our 
asynchronous setting. As a result, we compare its performance 
with existing methods in estimating the covariance matrix 
and AoAs. To test the MSE of our proposed covariance 
estimator S&M, we compare its performance with the sample 
covariance, shrinkage method in 1,24] and shrinkage MMSE 
||23]| estimates. For these tests, we use the temporal covariance 
matrix with d = M = 72. For all simulations, we set 
SNR dB, \Ci\ =6 for all i G I and let K range from 
5 to 50. Fig. [3] shows the normalized MSE of the estimators 
defined as 

lis -Sill 

mil ' 

for different values of K, where S denotes any of the covari- 
ance matrix estimates. It is evident that S&M outperforms all 
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Fig. 2. pT and pj versus SNR when Nr = 12 for K = 20 and K = 30. Fig. 4. AoA estimation RMSE when K = 20 for Nr = 10, 12, 14. 



the methods and in addition, it is very robust when K varies 
and even for iiT = 20 an acceptable performance has been 
achieved. 
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Fig. 3. MSE of covaiiance estimators for 100 channel realizations. 

In Fig. m the performance of S&M in estimating the AoA 
is presented in terms of root mean square error (RMSE) and 
compared with root-MUSIC (Muhiple Signal Classification) 
algorithrrQ. The curves are obtained by averaging the results 
of 1000 independent simulation runs for different SNR values 
and N]^. The RMSE decreases by SNR and the gap between 
the root-MUSIC and S&M RMSE reduces by increasing the 
number of antennas. However, in low SNR values, S&M 
outperforms root-MUSIC even for large number of antennas. 

VII. Conclusions 

We proposed a new algorithm for spatial and spectral 
sensing in asynchronous MIMO-OFDM signals. Our method 

' Even though they are not shown here, we have compared the performance 
with standard MUSIC and ESPRIT. However, the root-MUSIC algorithm 
had the best performance for AoA estimation among the methods based on 
unstructured covariance matrix estimation. 



first finds an accurate estimate of the covariance matrix of 
the received signal at the SU using shrinkage method for 
small sample size. Then, it matches the estimate with a sparse 
representation of the covariance matrix which depends on the 
channel parameters and occupied subcamers. The SU finally 
uses the support of the estimated sparse coefficient vector to 
recover the spatial and spectral pattern of the received signal. 

Appendix A 
Proof of Lemma[T] 

In the asynchronous scenario, at most two consecutive 
OFDM symbols transmitted from the ith PU are captured in 
each multipath contribution of y^, r = 0, . . . , K — 1. The 
NbM X 1 vector of received signal for any of the K collected 
observations is written as = + w^, where is the 
received signal due to the transmitted signals from the PUs and 
Wr is the noise. The vector of received signal corresponding 
to the transmitted signals x^ can be written as follows 

L 

X. =5]^x5;(z,£)+xi(z,€), (22) 

where at the rth observation period the contribution from 
the £-th path of the i-th user Xr(i,£) is written as the sum 
of the contributions of two consecutively transmitted OFDM 
symbols. The vector xj(i, l) corresponds to the OFDM symbol 
which its beginning is captured and x"(i,^) denotes the 
previous OFDM symbol which its tail is being observed. 

Each of the vectors x.l{i,£), j = 0, 1 can be written in 
terms of the channel parameters and the random transmitted 
data symbols as 

(23) 

where the indices kj, j = 0, 1 are the indices of the observed 
OFDM symbols at the rth observation sequence corresponding 
to the ith multipath component of ith PU. They can be 
described as follows where for simplicity we have dropped 
their dependency on i, £ and r in their symbol names 

fco = 2r - C - 1, fci = 2r - C, (24) 
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where 



M 



The parameter Vi^^ <E A4 is, defined 



as Vi^e ^ {ti + qi,i)[mod M) and represents the relative 
displacement between the ^-th multipath component of the 
observed OFDM symbols transmitted from i-th user and the 
start of the observation window at the SU (see Fig. [T]!. The 
vector aj^' = vec(A*^' ), j = 0, 1 is the vector of data symbols 
of the ith source in the frequency domain where non-zero rows 
in Aj^ correspond to the active subcarriers of the ith source. 
The matrix J^ (wi^^) is the shift matrix and is defined as 



3\v) = 



(25) 



Ot,xM-t) 

Odx Af-D 

Ia/— u 

The covariance matrix of the received signal at the SU is 
defined as S = E{y,.yf }. Since, the transmitted data and 
noise are independent, we can write 

S = Sx + I, 

where = EjxrX^}. In addition, a'i" and are 

independent for every i and by exploiting the fact that 

^{hliih^', I,,)*} = ai,eS{i-i')S{e-£')S{k-k'), we can write 
Ex = + where 

L 

-EEl'^W(*'^)(^'(*'^))''}- (26) 

iei i=i 

For notational convenience, we introduce the coefficient 

/ L 

"L,e = E E '^u.i ^e,e^Av - VuA^[i - (27) 
ti=i 1=1 

where v <^ M. and Q G A@ to indicate whether the ith source 
is being observed and whether there exists a link at a certain 
delay plus misalignment vT with a certain channel parameter 
B. The indicator 5c,c is defined as 

r A I li C 

Oc.c. = 



Otherwise ^'^^^ 
Using the parameterization introduced in (|27t , it follows that 

^Ahv.Q) = . (J-'"(z;) ® l)(r({*(m,0)}^froi)) 



(29) 



We define Oi^^^ = ^{^\ v e^'-^\ v e)*^ '^^e variance of 
the channel coefficient which in our setting does not depend 
on 2 and remains the same for all K collected observations. 
The expression for taking to account the parametrization 
introduced in ( |27] l is rewritten as follows 

" E E ^)(^' ^))''} 

Af-l 

= E E E '^^".^ ' (30) 

where 

0l)r({*(™,0)}ito')(F^diag({7.,c}lV)F* ® I) 
xr({*^K0)}it„l)((J■'■(^;))^®I), (31) 



where we have replaced E|aJ^^ (a*^^ )^| = diag({7i^c}^o^)® 
Iat^ in the expression. The coefficients 7^ c models whether 
the subcarrier c is occupied or not and is defined as 

'^^^ (32) 



otherwise 
The term F'^diag({7,;,J^^^)F* is simpHfied as 

F^diag({7„e}l-o')F* = E ^'.'^ f^f-^ (33) 
where fc = [F-'^Jc+i. Then, we can rewrite ( |30t as 



Af-l JV-1 

- E E 



^-^EE E E '^^-.^ 7^,c n^,,e ^ , (34) 



where 11;^, e c redefined as 

nl.,e = (J^"M S5i)r({*(m,0)}*to')(f,ff ®i) 

X r({*«(m, 0)}^-o^)((J^' W)^ ® I). (35) 

The expression in (l35T l illustrates that we can define the 
dictionary components as vec(n" ^ ^ + e c) which do not 
depend on i. Thus, as far as the estimation of the covariance 
matrix is concerned, we can define the new coefficients 

<yv,e,c — E 7i,c , (36) 

by summing over i. Using the vectorization operator, the 
sparse representation of the covariance matrix is more com- 
pactly written as 

vec(S) = Mo- + cr^ vec(lA/jvH), (37) 

where M is the over-complete dictionary and its columns are 
constructed with vec(n%^ + g^) for v e M., e A& 
and c G C The vector cr is an MNn^ x 1 sparse coefficient 
vector with non-negative entries equal to ay^e,c- 

Appendix B 
Proof of Corollary[T] 

In this case, Oi^i = [pi.^, &i,£, Qi.^], the corresponding dis- 
cretized vector is denoted by = [p, 6, a] where p ^V,a ^ A, 
h e B wd A& = V X B X A. 

In order to explicitly describe the dictionary components, 
we replace *(m,0) in (O with e^^''^e,.(6/-B)ef (a/A). 
Then, it follows that 

r({*(m, e)}^ll) = Aip) ® erib/B)ef{a/A), (38) 

where A{p) = diag e"* '^^/"o ^ is a diagonal matrix 
containing the Doppler coefficients. 

substituting ( l38b in ( |35] ), results in the following expression 

for K,pAa.c 

KpAa,c = (J^'(t') ® I)(A(p)fcff A^b) ® e,(VS)ef (aM) 

xet(aM)ef(&/S))((J-'(t;))^®I) 

= 7Vt(J^ WA(p)fcff A^(p)(J^(«))^) 

(g)er{b/B)e^{b/B), (39) 
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which does not depend on a. 

We define the matrix ^ ^ ^ which does not depend on the 
parameter a anymore as 

3l.,P.c = NT[3\v)K{p)U^K"{p){y^{v))"). (40) 

By the entrywise computation of (|40] l. one can verify that the 
sum H"^,.p_c + '^4>,p,c is simpHfied to 

+ Ht,p,e - Nt[Hv%^c ^" {V)) (41) 

where T(t;) denotes a masking matrix and is defined in d^TT l. 
The covariance matrix is eventually expressed as 

= XI H cr^,v,p,b %,c (*p,c r{v)) 

i^X v,p,b.c 

®{er{b/B)e^(b/B)), (42) 

where *p,c = A{p)iJ^ A" {p). The expression in (|42] | 
illustrates that we can define the dictionary elements as 
*p,c0f (^) in the temporal domain and er{h / B)e^ [h / B) in 
the spatial domain which do not depend on i. Similar to ( [36] l, 
we define the new coefficients 

crv,p,b,c = X Nt 0-i,v,p,b li,c- (43) 

In matrix form, the sparse representation of the covariance 
matrix is similar to ( iJTl l. The results of Corollary [T] follows 
where the vector cr is a MNPB x 1 sparse coefficient vector 
with non-negative entries defined as p ^,6) — <^v,p,b,c- 

Appendix C 

Derivation of the Spatial Covariance Matrix 

The N]^ X Nfj covariance matrix in the spatial domain is 
defined as Ss = E{y[m]y^[m]} = E{x[m]x^[m]} + a^I. 
We first find the covariance of the mth and the m'th received 
vectors at the SU as 

ieX 1=1 fc,fc'GZ 

xef(6,,,/i3)E{ef(a,,,M)A,fe4_,A,_,^_^_,^ 

X K"-k'M-,^,,-tA"k'et{aui/A)}. (44) 

The term {ai^i/A)Ai,ki^-kM-q- t-ti ii^ '^^e expectation in 
(|44] | is a scalar. As a result, one can rewrite the expectation 
as 

K"-k'M-g„,-U^{^"k'et{a.,,i/A)ef{a,,,/A)A,,k} 

X tm-kM-qi i-ti- 

By exploiting the independence of the rows of we have 

E{Af^k,et{a,j/A)ef{a,j/A)A,^k} 

= Nt diag({7,- - k'], (45) 



where 7,; c is defined in (l32T l. Thus, one can rewrite (l44l l as 

E{x[m]x^[m']} '^^-^ e^'2''''''%7"" 

iGl f=i fcez 

X f^-feA/-gi,f-tidiag({7i^c}^o^)fm-*:M-gi_f-ti 

X erib,j/B)e^{h,j/B). (46) 

The expression in (|46] | demonstrates that the mth and the m'th 
observations are uncorrelated if m' = m + kM for fc G Z and 
thus to collect K i.i.d. samples one has to choose m — kM 
for k ~ 0, . . . , X — 1. Specifically, when m — m', we can 
rewrite (l46t as 

E{x[m]x^H} = EE '^-^^^(l^)^^(^)- (47) 

In order to express the covariance matrix over an overcom- 
plete basis, we introduce the coefficient 

I L 

^^'b = E E ''^'^ - ^J.^l^t^ - -^l- (48) 
i=i fci 

Then, we can express S5 equivalently as 

Ss = E e,(VS)ef (VS) + I, (49) 

h=0 

where ^6 = matrix form, the sparse representa- 

tion of the spatial covariance matrix is more compactly written 

as 

vec(Ss) = Mso-s + (T^vec(I), (50) 
where [Ms],+i = vec(e,(VB)ef and [cTs]b+i = 6- 

Appendix D 

Derivation of the Temporal Covariance Matrix 

The M X M temporal covariance matrix for the ^th set 
of spatially filtered observations Yt(j') = €>^YT(r) in the 
temporal domain is expressed as St.; — E{y^;y*;} where 
Yr.i is the Zth row of Yrir). In the asynchronous scenario, 
at most two consecutive OFDM symbols transmitted from the 
zth source are captured in y,,; = (piYxir) where 4>i is the ^th 
row of and Yrir) = X,. + W^. The Nr x M block of 
received signal from the sources for any of the K collected 
observations can be written similar to (l22t as follows 

L 

iei e=i 

where vec{Xl{i, ^)) = x^(j, i) for j = 0,1. 

For notational convenience, similar to (l27t . we introduce 
the coefficient 

/ L 

^iv,p,a,b = E E '^nU - auMb - buA 
u=l £=1 

X S[v - Vu,e]S[p - puj]S[i - u], (52) 
where v G M, a e A and b e B. 
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Using the parametrization introduced in (|52t . one can ex- 
press X:^(i, £) as 

Xiii, v,p, a, b) = al^ p „ f,er{b/B)ef{a/A)S^i, v,p), 

(53) 

where S^{i,v,p) is defined as follows 

S\t,v,p) = A^-^FA(p)(P(t;))^. (54) 

The temporal covariance matrix St,; — ^{yJiyri) — 
E{X^0f 0*X*} +0-^110^1111 can be equivalently expressed 
as St,; =S?,^z + i:T,i+o2||0f||2l where 

■i^X v.p.a.b 

= E E TTb,; 7j,c ^i^v,p,c ' (55) 

and (Tiy^ab — lEja'' (,(q!'' is the variance of 

the channel coefficient which does not depend on j and is the 
same for all K collected observations. The matrix aj^ ^ ^ is 
defined in ( |40] |. In ( fSSl ). we have replaced 

E{ (A^)^ (aM)e^(6/i?)0r0r e: (VS)ef(aM) } 
= iVT ^6,; diag({7,, J^^-qI) , (56) 

where E{(A^)^e*(a/A)enaM)A^^} = diag({7.,c}t~o') 
and the parameter tti,j = \e^{b/B)4>i \^ is a scalar adding 
more sparsity to our model, because it is equal to zero when 

h e A4,. 

The expression in (|40] | illustrates that the dictionary ele- 
ments, (Sg p 2 + p ^), expressed in (l4Tl i. do not depend on i 
and a. Thus, as far as the estimation of the temporal covariance 
matrix and the occupied subcarriers are concerned, we can 
remove a from the indices and define the new coefficients 

B-l 

'^v,p,c — E E '^i'^'P''' li,c- (57) 

Using ( [57] l and (1411 1. the temporal covariance matrix is ex- 
pressed as 

St,; = E<p,c (*p,.0T(«))+a2||0nill- (58) 

v.p.c 

If hi G A,f, is the lt\\ estimated AoA with its corresponding 
steering vector as the /th column in $ and hi = hi^i, then the 
coefficients cr[ are nonzero for v — Vij, p = pij, c £ 
d. Thus, the filter $ makes the representation even sparser 
and more structured. The sparse representation of the temporal 
covariance matrix is more compactly written as 

vec(ST,i) ^Mto-t,; +(72 ||(/,f||2vec(lM), (59) 

where [o'T,i\p(v,p,c) = '^vpc is a MNP x 1 sparse coef- 
ficient vector with non-negative entries and \^t\p{v,p,c) — 
vec(*p,c0T(w)). 
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